/* 
 * File:   main.cpp
 * Author: jonk
 *
 * Created on 21 de noviembre de 2011, 14:06
 */

#include <cstdlib>
#include <complex>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>
#include<iomanip>
#include <sstream>

//#define CALC_STRESS
//#define POLYDISPERSE
//#define CONFINED
//#define CORRELADORES
//#define PRINTXYZ
//#define FLUJOS

using namespace std;

#define kB 1.0
#define pi 3.1415926535898
#define avogadro 602.205

using namespace std;

template<typename T>inline T SQR(T x)
	{  return(x*x);  }
template<typename T>inline T POW3(T x)
	{  return(x*x*x);  }
template<typename T>inline int anint(T x)
	{  return(int(floor(x+0.5)));  }

#include "utils.h"
#include "particulas.h"
#include "utilsPart.h"
//#include "caja.h"
//#include "correlador.h"
#include "interaccion.h"
#include "MD.h"
//#include "verletClass.h"
//#include "verletBrownian.h"
//#include "espreso.h"
//#include "print.h"



/*
 * 
 */
int main(int argc, char** argv) {
   ///Particulas
    particula *part;
	///Tiempo de la simulacion
	double t;
	///Parametro gamma de la resolucion de las ecuaciones de movimiento
	double gamma=10.0;
	///Energia potencial
	double epoten;
	///Empleo del termostato
	bool termo = false;
	///Caja de la simulacion
	//caja *box;
	///Correlador
	//correlador *correlator;
	///Cargamos el archivo de configuracion
	FILE *fe;
    /// Abrir el fichero de entrada en lectura y binario
    fe = fopen("config.bin", "wb");
    fclose(fe);
	parametros param;
	importConfig(&param);
	printConfig(param);
	///Nombre del archivo de tensor de esfuerzo
	char name[50];

    ///Semilla aleatoria
	//srand (time(NULL));

    ///Creamos las particulas
    part = new particula[param.npart];

	///Creamos la caja de simulacion
	//box = new caja(param);
	setBox(part, &param);

	if(param.rest==0)
	{
	    ///Comenzamos la simulacion desde el principio
	    SetRandomPart(part, &param);
	    setVelocidades(part, param);
	}
	else
	{
	    ///Reiniciamos la simulacion
	    loadConfig(part, param);

	}


	///Creamos el archivo binario de almacenamiento
	ofstream ffout("posvel_bin.dat",ios::out|ios::binary);
//	ffout.write((char *)(&box->npart),sizeof(int));
//	ffout.write((char *)(&box->lbox[0]),sizeof(double));
//	ffout.write((char *)(&box->lbox[1]),sizeof(double));
//	ffout.write((char *)(&box->lbox[2]),sizeof(double));
//	ffout.close();
//
//    #ifdef CORRELADORES
//	///Creamos los correladores
//	correlator = new correlador(box,param);
//	///Inicializamos el fichero de salida de los correladores
//	correlator->initFileCorr();
//	#endif
//
//	//verlet *vrl= new verlet(box,Dt,termo);
//	//verlet *vrl= new verlet(box, param.Dt, gamma, termo);
//	double temp=1.0;
//	//verlet *vrl = new verlet(box, param.Dt, temp, 0.05);
//	//verletBrownian *vrl = new verletBrownian(box, param.Dt, temp, 0.005);
//	verletEspresso *vrl = new verletEspresso(box, 10.0, 2.0, param.Dt);
//
//    #ifdef CONFINED
//    ///Calculamos la interaccion del sistema con confinamiento
//	interaccionConfined(box);
//	#endif
//
//	#ifndef CONFINED
//	///Calculamos la interaccion del sistema sin confinamiento
//	interaccionSimple(box);
//	#endif
//
//    ///Creamos el archivo de energia
//    ofstream fout ("energia.dat",ios::out);
//	fout.precision(5);
//	fout.is_open();
//
//
//	///Comienza el bucle principal
	t=0.0;
	for (int i=0;i<100000000;i++)
	{
//
		///Actualizamos el paso de tiempo
		t+=param.Dt;
		param.t=t;
		//cout<<"Ini  "<<part[1].vel[0]<<"  "<<part[1].posOld[0]<<endl;
		epoten = MD(part, param);

//		box->nacc=i;
//		//vrl->PaulYoonEngine();
//		//vrl->verletLangevin();
//		vrl->integrate_vv();
//
//		#ifdef CORRELADORES
//		correlator->calcCorrelador(t);
//		#endif
//
		if(i%10==0)
		{
//
//			//time_t tAct = time(NULL);
			printf("%10.5f  %8.2f \n",t,epoten);
			 printBinary(part, param, t, i);
		}
//			fout<<t<<"  "<<box->epoten<<"  "<<box->ekin<<"  "<<box->epoten+box->ekin<<"  "<<box->temp<<endl;
//			//printConfigBackup(box,param);
//			box->printConfigBackup();
//
//			#ifdef PRINTXYZ
//			box->printXYZfile("system.xyz");
//			//printXYZfile(box);
//			#endif
//
//			#ifdef CALC_STRESS
//			box->printShearStressTensor();
//			//printShearStressTensor(box, param);
//			#endif
//		}
	}
//    fout.close();
//	delete []box;
//	delete []correlator;
	delete []part;


	cout<<"sacabo"<<endl;
	system("PAUSE");   ///Particulas
        return 0;
}

